def zProb(x, type = 0):
    """Calculates the area under the curve of the standard normal
    distribution. (A)

    INPUTS:
    z (float): z-statistic
    type {int, 0}: {0,1,2} (See (1))

    NOTES: 
    (1) 0 = area under the curve to the left
        1 = area under the curve to the right
        2 = two-tailed test

    REFERENCES:
    (A) Algorithm AS 66: Applied Statistics, Vol. 22(3), 1973
    """

    x0 = 0.398942280444; x1 = 0.39990348504 
    x2 = 5.75885480458; x3 = -29.8213557808
    x4 = 2.62433121679; x5 = 48.6959930692
    x6 = 5.92885724438
    x7 = 0.398942280385; x8 = -3.8052e-08
    x9 = 1.00000615302; x10 = 3.98064794e-04
    x11 = 1.98615381364; x12 = -0.151679116635
    x13 = 5.29330324926; x14 = 4.8385912808
    x15 = -15.1508972451; x16 = 0.742380924027
    x17 = 30.789933034; x18 = 3.99019417011
    lt = 7.
    utz = 18.66
    con = 1.28
    pos = x >= 0
    z = abs(x)
    if z > lt:
        pvalue = 0.0
    else:
        y = .5 * z * z
        if z <= con:
            pvalue = .5 - z * (x0 - x1 * y / 
                    (y + x2 + x3 / (y + x4 + x5 /
                        (y + x6))))
        else:
            pvalue = x7 * np.exp(-y) / (z + x8 + x9 /
                    (z + x10 + x11 / (z + x12 + x13 /
                        (z + x14 + x15 /
                            (z + x16 + x17 /
                                (z + x18))))))
    if type == 0:
        if pos:
            pvalue = 1. - pvalue
    elif type == 1:
        if not pos:
            pvalue = 1. - pvalue
    else:
        pvalue = 2.0 * pvalue

    return pvalue

print zProb(5.8586497, type = 2)